On an "interaction by moments" property of four center integrals. 
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Abstract 

The four center integrals needed in the Hartree Fock approximation and in TDDFT linear response are known to be 
difficult to calculate for orbitals of the Slater type or of finite range. We show that the interaction of pairs of products 
that do not mutually intersect may be replaced by the interaction of their moments, of which there are O(N). Only 
quadruplets of orbitals 'close' to one another need an explicit calculation and the total calculational effort therefore 
scales as O(N). We provide a new and concise proof of this "interaction by moments" property. 



Motivation 

This note is motivated by the occurrence of four center integrals in the TDDFT linear response equation [T|, [2] 
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where Xks ~ sv^s ant ^ ^ = $v P " are ^ ne ^ ee anc ^ interacting density response, where p denotes the electronic density 
and where V ex t is an external potential acting on the electrons. In a basis of local orbitals {/ a (r)}, the electronic 
density p and the external potential V ex t may be expanded in terms of products of such orbitals according to 

p = E rt Pab , v& = J drv^rf (2) 

a, b 

Equation fTJ) then turns into a matrix equation for the response % a b cd = tvek °f P a b(fi) with respect to variations of 
Vext(t') and this equation then contains the Coulomb interaction between products of orbitals or 'four center integrals' 

< 12|i|34 >= / drdr'fM/ 2 ^)^-/ 3 *^')/ 4 ^) (3) 
r J \r — r | 

A technique for calculating these quantities in terms of two center integrals was developed in [3]. In an alternative 
"resolution of identity" method, products of orbitals are replaced by auxiliary functions, see [4] , [5] . 

For orbitals of finite range, there are, for N atoms, 0(N 2 ) pairs (12), (34) of individually intersecting orbitals 
suggesting the need of 0(N 2 ) distinct calculations. Here we show that only the subset of < 12|i|34 > where a pair 
(1,2) intersects with a pair (3,4) must be calculated explicitly, while the remaining ones can be taken into account 
by their multipolar interactions. Since there are only O(N) such quadruplets of orbitals and because the effort of 
calculating the moments scales like N, the cost of calculating four center integrals then scales as O(N). 

The present note arose in an ongoing effort to implement linear response for extended molecular systems, an 
effort prompted by recent work 0. We prove and exploit an observation of Greengard on the exact character of the 
"interaction by moments". This observation was rederived previously in the literature [6] and its consequence has 
recently been absorbed in a corresponding computer code [8], but our concise and simple deduction of this important 
property may still be of interest. 



Reduction from four centers to two centers 

We first need a reduction of products of orbitals to a set of single center functions. Following the discussion of [3] we 
obtain an expansion in spherical harmonics of a translated function /( r — a ) by using its momentum representation 

A m C?-^) = / J0j2^ m (lt)e- ip(r - a) (4) 
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Figure 1: convergence towards a cusp with j max 

We use spherical Bessel functions ji(x) with I integer, that are related to conventional Bessel functions for half integers 
by jl(x) — \/^Ji+l( x )- ln practice, fast Hankel transform routines [9] are needed to speed up the calculation. 

Expanding e 1 v ' r in spherical waves, one finds 

">Pim(~r ~ ~a ) = Yi i' n i (^) G ii'"i ( r ) ( 5 ) 

where the coefficients G; imi (r) that multiply the spherical harmonics Yi imi ( r ) now depend on both quantum numbers 
h and mi because spherical symmetry is lost. One has the following expression for G; imi (r) 

G hmi (r) = 2v^^F„ li2 (r)y4 mi _ m (^)(-)' 1 (-)^ i2+il+i) G im , iimii2mi _ m (6) 
h 
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where Gi imi ,i2m 2 i3m 3 are Gaunt coefficients for the overlap of three spherical harmonics. To bring this approach into 
perspective, it is interesting to consider an orbital with a cusp singularity such as e~ r and to study the convergence of 
the translated orbital towards a translated cusp with increasing angular momentum cutoff j max , see the figure. The 
figure shows that a fairly large number of angular harmonics is required for the representation of orbitals having such 
a cusp. 

By applying translations to two overlapping orbitals we may obtain an expansion of their product about a common 
midpoint 

V' !imi (^ > -"ai)V'( 2m2 (^"""2) = Gim(r)Yi m (~r) (7) 

i=0..i max , m=-l..l 

Using such expressions we can compute the Coulomb interaction of two pairs of mutually overlapping orbitals from 
the point of view of their associated effective centers. 

Exact interaction by moments of non intersecting centers 

We have seen that pairs of orbitals may be replaced by effective centers. Now we wish to show that the interaction 
between effective centers may be simplified when their spheres of support no longer intersect. We begin by quoting 
formulas for the computation of two center integrals in Fourier space, see [3], [10] 




where we used that \ ri l T2 \ = A being the Laplace operator. Expanding e 1 p ' r in spherical waves one finds 
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with 



C lir iMm ^{ R ) = Cl x i 2 i(R)Gi imi ,l 2 m 2 ,tni (9) 



where G; imii ; 2m2 ,; m are the previously encountered Gaunt coefficients and C; 1 ; 2 ;(/?) are Wigner-Eckart like couplings, 
see [3], [10] for details of the derivation. We use eq(|4| to rewrite the Ci lt i 2t i(R) in terms of the original radial wave 
functions as follows: 

l+h^h f 2 ~» 2 ~ f°° 

C lu i 2 ,i = 32ir(-) 2 J dr 1 dr 2 r 1 i[: ll (r 1 )r 2 ip l2 (r) J ji 1 {pri)ji 2 {pr 2 )ji(pR)dp (10) 
The "interaction by moments" property we are after is contained in the following integral of Bessel functions 

/■oo 

L i a ,(ri,r 2 ,R) = / jh(pri)ji 2 (pr 2 )ji(pR)dp (11) 
./o 

The Coulomb interaction coefficients Ci lt i 2] i would reduce to I ; , 2 l (ri,r 2 ,R) if the original orbitals functions were 
concentrated at, respectively, radii n and r2. This integral therefore represents the Coulomb interaction of two 
charged hollow shells of radii ri , r 2 at a distance of _R, and with the charge densities having the appropriate multipolar 
angular dependences. For R > r± + r 2 where these shells no longer intersect, we expect their interaction to simplify. 
Because of the Gaunt coefficients in eq© we only need this interaction for even values of h + h + I and where a 
triangle inequality \h — l 2 \ < I < h + l 2 holds. In this case the integrand associated with the dp integration in eq (|10p 
is symmetric as a function of p and we may therefore extend the domain of integration to the entire p axis: 

1 f°° 

I i 1 ,i 2 ,i( r iy r 2,R) = y / jh(rip)ji 2 (r 2 p)ji(Rp)dp for h + l 2 + I even 
J —00 

It is convenient to use ji(z) = Kehi(z) and to consider a corresponding complex integral / ; c ; ( with I f ia ; = Re/ ; c , ; 
that involves the Hankel function hi(z) : 

1 f°° 

I^ l2l (ri,r 2 ,R)= -= \ ji 1 {r 1 p)ji 2 {r 2 p)hi{Rp)dp, (12) 

For R > n + r 2 the contour of integration in 1° ^ l (n, r 2 , R) can be be closed at infinity in view of the relation 

p 

(13) 



Clearly, the exponential factor e tHp hom ji(Rp) dominates, for R > ri + r 2 , the factors e ±lpri e ±lpr2 that arise in the 
integrand in eq (|12p from the product jn 1 {rip)jn 2 {r 2 p). Since the integrand in 7 ; c ; ; is analytic, except for possible 
singularities at p = and since the contour of integration can be closed in the upper half plane, a non zero contribution 
to 7° t , can om y be due to a residue at p = 0. From eq (|13[l the most singular term in hi(Rp) at p — is 

hl (Rp)=-i(2l-l)l\^^+0(p- 1 ) 

Because of j^irip) ~ (2^+1)!! an< ^ an analogous relation for ji 2 (r 2 p) a non zero residue is impossible in eq (!12p unless 
I attains the maximal value permitted by the triangle inequality, I = Zi + l 2 . When setting I = ii + l 2 and closing the 
contour of integration in eq (|12[) at infinity, there is a term ~ 1/p that provides a non zero result by elementary contour 
integration. Rewriting the result in terms of conventional F functions, one then finds, for I + l\ + l 2 even and with 
< I < h + h, the following simple result 

( m = h 7r 3/2 r^r^ £([+1/2) 

h^A^^R) d wi+i2 8 Rl+1 r(/ 1+ 3/2)r(; 2 + 3/2) 1 j 

Applied to the Coulomb interaction coefficients of eo (|10p one concludes 

w^ii tv 1 ' 2 P ' h P l2 £([+l/2) 
^ 2 ; i+ ; 2-'+2 ff+i r(Z 2 +3/2)r(Z!+3/2) Ul+i2 1 ; 

fti 2 = 4?r / drr 2+ll - 2 i) ll (r) iorR>n+r 2 
Jo 

This last equation shows very clearly that non overlapping orbitals interact exactly via their moments, as shown first 
by Greengard [11]. For another proof, see [6]. 
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Conclusion 



We conclude that only the subset of four center integrals < 12| 1 |34 > of "close" pairs where (1, 2) intersects (3, 4) must 
be calculated explicitly. Because there are only O(N) such pairs of orbitals for a system of N atoms and because the 
multipoles are associated with only O(N) products, the calculational effort scales as O(N). 

The conclusion that Coulomb integrals should be divided into far and near field ones has already been incorporated 
in a quantum chemistry code [8]. But our derivation of the "interaction by moments property" of four center integrals 
from a plain integral of a product of three spherical Bessel functions is the simplest and most concise proof of this 
property that is available. 
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Figure caption: Translation away from the origin of e~ r by one unit, with j < j max = 4,6,8,10 and the 
convergence of the result towards a cusp with increasing angular momentum cutoff j max . 
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